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Abstract 

We implement a Quantum Monte Carlo calculation for a repulsive Hubbard 
model with nearest and next-nearest neighbor hopping interactions on clus- 
ters up to 12 X 12. A parameter region where the Fermi level lies close to 
the van Hove singularity at the Saddle Points in the bulk band structure is 
investigated. A pairing tendency in the d^2_y2 symmetry channel, but no 
other channel, is found. Estimates of the effective pairing interaction show 
that it is close to the value required for a 40 K superconductor. Finite-size 
scaling compares with the attractive Hubbard model. 
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The 2D Hubbard model (HM) contains several basic elements of the high temperature superconductivity 
problem, and its properties include nontrivial features (e.g. the Mott transition and associated antiferro- 
magnetism) which are generic in HiTc materials. Is the HM also a superconductor, or does the supercon- 
ductivity originate from extrinsic interactions or degrees of freedom? According to calculations , |^ based 
on exchange of antiferromagnetic spin fluctuations (AFSF), the HM is a superconductor with d,j.2_y2 order 
parameter symmetry (a symmetry consistent with several recent experimental measurements |^). But 
on the contrary Quantum Monte Carlo (QMC) calculations up to the present have given negative or 

inconclusive results regarding superconductivity of the 2D HM. 

In the present QMC study of a generalized HM (which includes next-nearest neighbor hopping interac- 
tions as well as the usual nearest- neighbor ones), we present evidence for a superconducting tendency in the 
d,j.2_y2 channel in large (up to 12 x 12) clusters, in qualitative agreement with the the AFSF calculations. 
The results have been obtained close to the parameter region where the energy of the saddle points (SP) in 
the band structure lie close to the 'Fermi level' in the cluster energy level structure, a choice motivated by 
the hypothesis - 'Van Hove Scenario' - that, in the continuum limit, superconductivity is enhanced 

by the Van Hove density of states peak associated with the saddle points. 

The SP feature may be incorporated into the Hubbard model within the metallic regime by introducing 
a next-nearest neighbor interaction. This allows the SP to lie at the Fermi level at a doping of, say 15- 
25%, while the insulating point, at which the antiferromagnetic instability occurs, lies at 0% doping. These 
features are characteristic of real cuprate materials |]l3|, |l4|. It is in the former situation (15-25% doping) 
that the model is found to support superconductivity. 

The model is specified as follows 



In (1), U is the repulsive on-site Coulomb interaction, t is the nearest-neighbor hopping integral ((ij) denotes 
nearest neighbor interactions), t' is the next-nearest neighbor hopping integral {{{ij)) denotes next-nearest 
neighbor interactions); t and t' are defined to be positive. 
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The noninteracting band structure of the ti'-Hubbard model (1) has saddle pomts at energy — 4i', and 
at fc = (0, tt) and (7r,0). If we take the hole doping x to be of order x k, t' , then the saddle points in the 
noninteracting band structure lie near Ep (see Fig. 1). In the absence of t' , the required doping would be 
zero, making the sample insulating. The electronic effective mass below the SP's is heavier than the mass 
above, the ratio being (t + 2t')/(t — 2t'). Photoemission data indeed show a large mass ratio [|l3|, jl^. The 
simulation of the tt'-Hubbard model |^ is based on the projector Monte Carlo technique, using the Ansatz 
for the ground state 0] 

1*3) = e®^|*o);e^oo (2) 

where is a projection parameter and |^o) a single determinant taken as the ground state of the noninter- 
acting band structure of |l|. 

The exponent in (2) is broken up into Trotter slices on the interval 0, and the Hirsch-Hubbard- 
Stratonovich [^5| transformation applied to each slice, turning (2) into a path integral which is solved 
on a finite cluster (up to 12 x 12) using the Metropolis algorithm (for details see Ref. §). We checked 
on the convergence of the algorithm by comparing small clusters to exact diagonalization and stochastic 
diagonalization |^ and found = 8 and r — 0.125 to be sufficient. Periodic boundary conditions and 
closed shell configurations are always used. The 'average sign' is usually sufficiently close to unity, that the 
simulations are not significantly limited by the fermionic sign problem. In a typical simulation we averaged 
over four runs with different seeds each run having at least lO^MCS. The nonexistence of a fermionic sign 
problem in the simulations and the very large number of MCS allowed us to obtain precise results for the 
superconducting correlation functions. The orthogonalization technique originally proposed by Sorella et 
al. helped us to stabilize the algorithm. We found a stabilization every 8 Trotter-slices to be sufficient 
to make the results independent of the number of stabilizations. 

The tendency towards superconductivity is signalled ||l^ by the presence of a long-range 'plateau' 
in the superconducting correlation 'function' plotted as a function of distance. Such plateaus have been 
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demonstrated [[7| in models such as the attractive Hubbard model and simplified models with electron- 
phonon coupling, where a pairing tendency is anticipated to occur. While studying finite systems one has 
always to extrapolate to the thermodynamic limit. Due to the unsystematic finite level structure of the finite 
Hubbard model this has turned out to be very difficult. We use two different approaches: a) we are able to 
analyse the plateau in terms of the effective pairing interaction J, and hence deduce a value for in the 
infinite-sample limit, b) we find analogous finite size scaling behaviour for the superconducting correlations 
for the repulsive and the attractive Hubbard model. 

The condition as to whether the finite cluster of size L is 'superconducting' (correlation length ^ > L), 
or whether superconducitivity is suppressed by finite size effects < L), does not enter into these consid- 
erations and is not relevant for this paper. 

Some constraints on U are a) according to exact diagonalization calculations the Hubbard model is an 
insulator at half-filling, as observed experimentally, if U is greater than about lOt' b) the value U=2t 
has been used in calculations [|l8| of the quasiparticle lifetime broadening for the tt'-Hubbard model in 
the Renormalized Propagator Approximation, with results consistent with experimental data c) attempts 
to derive the single-band U formally from a multi-band model give a value of order 6eV, i.e. 6t, in 
the case t = 1 eV. We have observed clear evidence for the superconducting tendency for U in the range 
0.5t <U < 3t; for larger U-values unnacceptable error bars are obtained. 

The superconducting correlation function x{Rj) for di^2_j^2 symmetry is defined as 




i (P)M) 



(3) 
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where 
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Here, (p) implies a sum over nearest neighbors of 0, p = (±1, 0) and (0, ±1), for which ap is respectively 1 
and -1, Nl is the number of atoms in the cluster, and a is spin. 

A typical calculation of the superconducting correlation function (3) for a (12 x 12) cluster is illustrated 
in Fig. 2a. It is seen that, as a function of increasing distance, the correlation function in this symmetry 
approaches a 'plateau'. It is only in d2,2_y2 symmetry that we find this behavior. The error bars in the 
points, calculated as the deviation between runs with different seeds, are smaller than their diameter. 

A study of the variation of x^' with U is shown in Fig. 2b, based on 8 x 8 clusters, x^'' is seen to increase 
with U, and continues to increase steadily up to the largest value of U (U=3t) for which we have adequate 
statistics. 

The main problem with scaling our results is the influence of finite shell structure on the observables and 
the unsystematic behaviour of these shells with system size. Therefore the application of finite size scaling 
especially for weak U is extremely complex, as has already been observed for the case of the attractive 
Hubbard model |2|. 

In Fig. 3 we illustrate a plot of the value of x{R) averaged over the plateau region, X'^\ versus t' for a 
10 X 10 cluster with N,, — 74 electrons. Each point in Fig. 3 represents a run with several different seeds, 
typically 4. Fig. 3 contains points calculated using two different sampling algorithms, which are seen to 
agree quite well. The strong variation of x^' with t' seen in Fig. 3 is systematic, and is due to the shell 
structure of the cluster. The peaks in x^' are found to be associated with the near-degeneracy of several 
filled and empty cluster levels, i.e. a high 'density of states at Ep\ The dips in x^' coincide with regions 
of maximum gap in the single particle spectrum of the cluster between the highest occupied and lowest 
unoccupied levels. 

A more quantitative analysis of the effect of shell structure on pairing is possible by introducing a simple 
low energy pairing Hamiltonian [ pO[ 




(4) 
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Here k — a discrete wavevector for the cluster (periodic boundary conditions assumed) lying within the 
cutoff |£fe — mI < ^c, Pfe = c^^c^k^, Vk ~ {cos{kx) — cos{ky)), and J, which is negative, is a nearest-neighbor 
attractive interaction which parametrizes the pairing strength in the model (4). 

For the present purpose we exploit the cutoff — /i| < lUc, to solve (4) by exact diagonalization for states 
within ojc- Having available the exact system ground state we can calculate x^', which is also illustrated 
in Fig. 3 for J=0.055t. ft is seen that the low energy Hamiltonian (4), with its shell-structure effects, well 
reproduces the corresponding structure in the QMC results. Using the complete eigenvalue and eigenvector 
spectrum of the of Hp, within the cutoff, the projector formula for x^', Eq. (2), can be evaluated explicitly 
for a given 8; we find that the value Q = l/8t adequately reproduces the ground state value of x^"'. 

In Table I we demonstrate a possible technique for scaling QMC x^' data on different cluster sizes. The 
problem is to avoid the pronounced cluster shell-structure effects (Fig. 3). The technique is to compare 
J, obtained by fitting the Eq. (4) model to x^' at each cluster size, rather than for different sizes, 
whereby the shell-structure effects are taken into account. Values of J deduced from QMC calculations on 
n X n clusters with n from 6 to 12 are given in Table I. It is seen that the results appear to be converging 
on J=-0.15t, though results on a wider dynamic range of cluster sizes are needed to be convincing that 
convergence on the bulk limit has been achieved. 

The foregoing method of analysis is extensible to include the case of a constant Z-factor (Eliashberg no- 
tation). The analysis yields J/Z^ . A constant DOS po renormalizes to ZpQ, giving Z /{J p^) in the exponent 
of the the Tc equation, which is the correct form and hence does not require any additional correction for Z. 
It is therefore valid, to within the constant-Z approximation, to proceed as if in the BCS case (Z=l). The 
largest value of J we have found is in the case t' — 0.22, U — (see inset Fig. 3), for which J = —0.2t. 
This value of J is close to that (J ~ — 0.22i) required to give a 40 K Tc, with t = 1 eV. 

In Figure 4 we present a second method to extrapolate the data. We found in agreement with ||2^ that 
scaling laws of the form 1/N are too simple and can lead to ambigous results. Instead of incorporating 
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corrections to scaling while applying the scaling laws |E7[ we directly compare the results on x'J for 
positive U with results on xf' for negative U. In both cases we use the same system sizes, fillings, and t' . 
Hence the finite size structure is identical. Additionally we use a value of U=-0.3 to match the size of the 
plateaus for the positive U {dx2_y2) case. Our analysis shows that the attractive and repulsive HM have 
the same scaling behaviour. As the attractive HM is superconducting we conclude from the above that the 
repulsive HM is also superconducting. 

Further light on the mechanism of pairing in the HM comes from an independent study of the Hubbard 
model based on the Roth decoupling procedure p2| . This study shows pairing of d^^-y^ symmetry, and 
exhibits a clear correlation between pairing and nearest-neighbor spin correlation. Therefore the source of 
the pairing is probably a nearest neighbor antiferromagnetic exchange interaction, generated by higher order 
diagrams (23| in U, starting at U'^. 

In conclusion, pairing correlations in d^^-y^ symmetry, but no other symmetry, are found for the clusters 
and parameters studied in this paper. We presented two different ways of extrapolating our data to the 
thermodynamic limit. The tt' HM phase diagram thus includes the Mott transition, marginal Fermi Liquid 
behavior, and superconductivity - all the key features present qualitatively in real cuprates. The pairing 
interaction found in the parameter space explored so far supports a Tc of order 40 K. Full understanding of 
high temperature superconductivity, including the highest Tc's, probably requires consideration of additional 
interactions, such as superexchange from the oxygen bands and mediation by phonons. 

The Monte Carlo simulation was carried out on the IBM SPl. Typically the simulation at a single 
parameter occupied 4 processors for 80 hours. 
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FIGURE CAPTIONS 



Fig. 1 Fermi surface of noninteracting tt'-Hubbard model t' = 0.3t, Ep = — 1.2t, hole doping = 0.27. Solid 
circles indicate saddle points. 

Fig. 2a Superconducting vertex correlation function x(^) (Eq. (3)), plotted versus distance R, 12 x 12 
lattice, 106 electrons (doping 0.264), 6 = t' = 0.286f, U = 2t. Insert shows correlation function 
on expanded scale, horizontal dashed line is average plateau value. Error bars are less than the width 
of the points. 

Fig. 2b QMC calculation of plateau vs. U, for 8 x 8 lattice, 50 electrons (doping 0.22), 6 = 8/t, t' = 0.22 

Fig. 3 Filled points QMC calculation of plateau value x^' of superconducting vertex correlation function 
(Eq. (3)) versus t' , for 10 x 10 lattice, 74 electrons (doping 0.26), = 8/t, U = 2t, error bar is 
average; open points fbcs calculation (Eq. (4)) with J=0.055t and cutoff omegac = 0.2t. 

Fig. 4 Superconducting vertex correlation x^'' in the onsite-s channel for attractive U = —0.3 versus su- 
perconducting vertex correlation x^'' in the dx^-ys-chamnel for the repulsive Hubbard model {U = 2). 
Filling and t' are the same for both values of U and are indicated in table I. 
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TABLES 



36 26 0.3 1 1.377E-03 0.12 t 

64 50 0.25 t 0.648E-03 0.15 t 

100 82 0.25 t 0.491E-03 0.15 t 

144 122 0.25 t 0.332E-03 0.15 t 



Table I Scaling t' = 0.22i, U = 2t 
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